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Abstract 

We describe numerical techniques used in our construction of a 4th order in time evolu- 
tion for the full Einstein equations, and assess the accuracy of some representative solutions. 
The scheme employs several novel geometric and numerical techniques, including a geometri- 
cally invariant coordinate gauge, which leads to a characteristic-transport formulation of the 
underlying hyperbolic system, combined with a "method of lines" evolution; a convolution 
spline for radial interpolation, regridding, differentiation and noise suppression; represen- 
tations using spin-weighted spherical harmonics; and a spectral preconditioner for solving 
a class of 1st order elliptic systems on S"^. Initial data for the evolution is unconstrained, 
subject only to a mild size condition. For sample initial data of "intermediate" strength 
(19% of the total mass in gravitational energy), the code is accurate to 1 part in 10^, until 
null time z — 55M when the coordinate condition breaks down. 
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1 Introduction 



The Einstein equations present difficulties of size and complexity somewhat greater than those 
normally encountered in scientific computation. A numerical simulation must confront a range 
of theoretical and numerical challenges, starting with the theoretical problem of finding a well- 
posed and geometrically natural reduction of the full system of Einstein equations. Numerical 
algorithms are then needed to control the various facets of the reduced system, efficiently, since 
very large data structures are inevitable when modelling fully S+l-dimensional spacetimes. 
Finally there are the twin theoretical and practical problems of understanding the nature of 
the solution represented by the data, and of certifying its reliability. 

In we presented a new coordinate formulation for the vacuum Einstein equations, based 
on a characteristic (null) coordinate and a quasi-spherical foliation Q . This null quasi-spherical 
(NQS) formulation is well-adapted to modelling spacetimes containing a single black hole, ex- 
tending from the black hole to null infinity. 

The purpose of this paper is to describe the numerical algorithms we have used in our imple- 
mentation of the NQS Einstein equations, and to present results of some accuracy tests of the 
code. Interactive access to the data sets described here, and many other simulations, is available 
online at http : / / gular . Canberra . edu . au/ relativity . html. More detailed discussions of the 
physical and geometric significance of the results of the code will be presented elsewhere. 

From the numerical programming viewpoint, the most significant features of the code are: 



1. a characteristic coordinate z ~ t — r (cf. [56|) plays the role of "time", with the numerical 
evolution in the direction of increasing z] 

2. the numerical grid is based on spherical polar coordinates on the z-level sets Mz\ 

3. S"^ dependencies are handled by a combination of spectral coefficients with respect to 
a basis of spin- weighted spherical harmonics (for spins 0,1,2, corresponding to scalar, 
vector and tensor harmonics); values of the field on a uniform grid in the spherical polar 
coordinates (1^,93); and Fourier coefficients of the values on the polar coordinate grid; 

4. a non- uniform radial grid along the outgoing null hypersurfaces Nz, which compactifies 
future null infinity and is adjusted dynamically so radial grid points approximately 
follow the inward null geodesies; 

5. reformulations of the hypersurface (radial) Einstein equations which enable the numerical 
modelling of the asymptotic expansions of fields near null infinity; 

6. 4th order Runge-Kutta time evolution; 

7. the characteristic transport (hypersurface, radial) equations are treated as a system of 
ordinary differential equations and integrated using an 8th order Runge-Kutta method 



8. an 8th order convolution spline is used for interpolation, for computing radial derivatives, 
to realign the fields with the dynamically varying radial grid, and for suppressing high 
frequency modes; 

9. a first order elliptic system on S"^ is solved at each radius and time step by the conjugate 
gradient method, accelerated by a geometrically natural spectral preconditioner. 
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A number of numerical consistency checks suggest that most quantities of interest which are 
calculated by the code (eg., the NSQ metric functions, the connection coefficients and the Weyl 
curvature spinors) have relative errors of about 0.001%, for simulations where the gravitational 
waves carry no more than 20% of of the total spacetime mass. Of course, greater accuracy is 
found for more nearly linear, weak field, simulations. The major limiting factor in determining 
the accuracy appears to be the spherical harmonic resolution, currently at L < 15. Although we 
have implemented routines which extend this to L < 31, a full implementation is not possible 
due to limitations in our present hardware. 

The algorithm was initially developed and tested on a 300MHz DEC Alpha with 512Mb 
memory, and presently runs on a 300MHz Sun Ultra 2 with 784Mb, with a typical (L=15) run 
taking between 2 and 4 days. Preliminary descriptions of these results are given in M, W^, 5^, 



40|. 



The Cauchy and characteristic initial value problems differ significantly in the nature of their 
appropriate initial conditions. The Cauchy problem initial data consists of the initial 3-metric 
and extrinsic curvature |54|, whereas the appropriate initial data for a characteristic initial 
surface is just the null metric. In the NQS case the null metric is 

ds]^^ = {rdd + l3^drf + {rsindd^ + p"^ drf, (1) 

parameterised by the angular shear vector [i = -^(/^^ — = Piz,r,i} tp). Unlike the Cauchy 
problem, the NQS Einstein equations Q do not impose any additional algebraic or differential 
constraints on /?, so the initial data P{z = 0) is an arbitrary function of the spherical polar 
coordinates (r, -iJc/p), except for a mild size constraint (|l9|). Heuristically /3 represents the in- 
going gravitational radiation of the spacetime; an interpretation which is consistent with the 
initial (z = 0) values of (5 being freely specifiable. Note that the geometric invariant onp/ Pnp 
m becomes -2S/3/(2 - div/3) in the NQS coordinates [|. 

The NQS geometric gauge provides a formulation of the Einstein equations as an explicit 
characteristic transport system, coupled to a time evolution equation. This type of structure is 
also found in characteristic formulations of other hyperbolic equations such as the wave equation, 
and it is well-known for the Einstein equations in other characteristic coordinate based gauges 



|28, 5C, |56|, 37]. Although there are existence results for characteristic initial value formulations 



of hyperbolic equations, these rely on reduction to the Cauchy problem |48, 66 1 rather than 
directly on the transport form. The only exception of which we are aware is the analysis of 
the linear wave equation in |2|, ^. It would be valuable to have theoretical existence results 
for systems of characteristic transport equations, which could justify the numerical formulation 
described here. 

The characteristic-based approach has recently been strongly advocated by Winicour and his 



coworkers [12, 25], who have developed codes for solving the scalar wave equation ]26], axially 



symmetric spacetimes ]24] and for the full Einstein equations Jll], based on the Bondi coordinate 



system [|2^. These works have been fundamental in establishing the feasibility of numerical 
formulations based on a characteristic coordinate and in motivating the present implementation. 
However, the stability analyses and experience of l p6| , p4[ are not directly applicable to our 
evolution procedure, since the formulations and numerical methods used have several significant 
differences. Like the Bondi parameterisation, the gauge conditions are directly implemented in 
the metric form; however, the NQS Einstein equations are considerably simpler than the Bondi 
coordinate form of the equations |5|, pi] ]. 

The treatment of angular derivatives is greatly simplified by the quasi-spherical condition, 
which encourages the use of spherical harmonic expansions. These in turn will simplify compar- 
isons with theoretical results on perturbations of the Schwarzschild black hole [E^, 5^, f^. 
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We note that the power of spectral methods in practical applications is well-known, in fields 
such as meteorology [17|, astrophysics |13], and fluid dynamics Finally, the combination of 
the characteristic-transport and method of lines techniques, considerably simplifies the use of 
high order algorithms such as RK4 for time evolution. 

The method of lines approach to evolution equations is common in fluid dynamics [^] , but 
has not previously been attempted in numerical relativity. The situation we consider here, 
of smooth variations in a single black hole geometry, is well-suited to the method, since high 
frequency modes are less likely to be of physical interest and thus may be treated by smoothing 
or artificial viscosity. 

The use of higher order methods, together with spherical harmonics and radial smoothing, 
leads to considerably more accurate results than would be possible with the 2nd order methods 
more commonly employed. However, this suggests that our code is restricted to very smooth 
spacetimes, and cannot reliably treat spacetimes with strong localised features (such as plan- 
ets, or gravitational shocks). Because we are concerned primarily with the gravitational wave 
perturbation modes of a single black hole, this does not present an important restriction, since 
the dominant modes are known to occur only for low angular momentum / - in particular, the 
/ = 2, 3, 4 modes are expected to carry practically all the radiated gravitational energy. 

The code models evolution in an exterior domain, and consequently boundary data must 
also be prescribed at an interior boundary surface r = tq. In general only two pieces of the 
boundary data can be freely specified, thereby fixing the null hypersurfaces A4 and the out- 
going radiation flux at the boundary The rest of the boundary data is constrained by 
boundary evolution equations (see (p^) , (p7|) ) . For simplicity, the version of the code reported 
here assumes fixed boundary conditions, corresponding to the past horizon of a Schwarzschild 
black hole. This assumption considerably simplifies the treatment of the inner boundary, and 
is completely consistent with the Einstein evolution. It follows that the code models the inter- 
action of gravitational waves with a single Schwarzschild-like black hole. Future versions of the 
code will incorporate dynamical inner boundary conditions. 

The paper is organised as follows. 

Section 2 gives the NQS form of the equations, and describes the structure of the equations 
and the formal steps in the solution algorithm. The geometric significance of the resulting 
equations is described in 

Section 3 describes the numerical techniques used, including the representation of spin- 
weighted spherical harmonics used to encode the angular variation of the various fields, and the 
high order convolution splines used for interpolation and differentiation in the radial direction. 

Two aspects of the treatment of spherical fields appear to be non-standard |42, 17|: the use 
of FFT's in both the (p and -d directions, based on the "torus" model of S*^ |3£]; and the use 
of spin-weighted spherical harmonics to handle, in a unified and frame-invariant manner, vector 
and higher rank tensor harmonics as well as invariant derivative operators. 

Section 4 describes the various stages in the evolution algorithm — solving the hypersurface 
equations out to null infinity T"*" ; reconstructing the metric from the connection variables (which 
includes the solution of a 1st order elliptic system on the 2-sphere at each radial grid position); 
and the evolution of the primary field /?. 

In section ^ we describe various techniques for estimating the accuracy of the code, test- 
ing both numerical and geometric properties of the numerical solution. Numerical convergence 
tests estimate the effects of separately increasing the resolution in the radial, time and angular 
directions. The two geometric tests described here demonstrate the consistency of the numer- 
ical metric by verifying the constraint equations for the Einstein tensor components GnmGnm 
and by testing the accuracy of the solution at infinity using the Trautman-Bondi mass 
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decay formula 14, 28]. These provide highly nontrivial tests of the consistency of the 

numerical solution. 

Other tests of the code, based on geometric properties of vacuum spacetimes, and compar- 
isons with known exact and approximate solutions of the Einstein equations, are envisaged for 
future work. 

2 Einstein equations and NQS metric functions 
2.1 Spacetime metric 

We consider spacetimes admitting global null-polar coordinates {z, r, (p) in which the metric 
takes the null quasi-spherical form 

ds%Qs = -2udz{dr + vdz) + (r dt? + dr + 7^ dzf + {r sm.'d dip + 0^ dr + -f^ dz)"^ , (2) 

where u > 0,v, and /? = P^d^) + 0^ cscdd^p, 7 = 7^9^ + 7^csci?9^ are the unknown NQS 
metric functions, to be determined by numerical solution of the Einstein equations. Note that 
we use , 9^ to denote equivalently the coordinate tangent vectors and the coordinate partial 
differential operators We may consider /3, 7 either as vector fields on 5^ or as spin 1 

quantities, defined by the complex combinations |23] 



The canonical example of a spacetime with metric in NQS form is Schwarzschild spacetime 
in Eddington-Finkelstein retarded coordinates |30| 

dslchw = -^dz {dr + i(l - 2M/r) dz) + r^{d^f^ + sin^ d d^"^) (4) 

with u = 1, V = — 2M/r), = 7"^ = and M = const. This includes Minkowski space 
M^'^ as the case M = and z = t — \x\. 

2.2 Edth 

Using the complex notation (ph, we have the canonical angular covariant derivative operator 



"edth" [pa, 44, 21 



acting on a spin s field 77, and its "conjugate" operator 6, 

All geometric angular derivative operators may be defined in terms of 8, 8. For example, the 
covariant directional derivative of a spin s field r] in the direction (3 is 

the divergence and curl (of a vector) are 

div/3 = g/3 + = + V2/3^ (7) 

curl/3 = i (8/3 - g^) = V2/3^ - Vi/32 ; (8) 
and the spherical Laplacian is 

At; = (gg + gg)r? . (9) 

Further properties of edth are described in section ^ and in 1 21 , 44 , ||] . 
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2.3 Connection variables 



In addition to the metric functions {u, v, /3, 7) we introduce the connection fields H, J, K, Q, 

H = -(2-div/?), (10) 
u 

J = v{2- div/3) + div7, (11) 

K = (12) 

Q = - + 7 + V;37 - V,/3, (13) 

g± = -{Q±du). (14) 
u 

Observe that u, v, H, J are real and have spin 0, whereas (3, 7, Q, Q^, have spin 1 and K has 
spin 2. 

Given the metric functions u, 7 on a z- level set Mz, we may construct H,J,K on Mz 
directly, and Q (and Q^) may be reconstructed if in addition, d/3/dz is also known on Mz- It is 
clear from (p!o|-[l^ that this construction does not require any compatibility conditions on the 
data u, V, /3, 7, g^p. 

Rather remarkably, there is a converse construction for the metric functions u,v,'j and 
df3/dz, which also involves totally free and unconstrained data, namely the connection variables 
H, J,K,Q. This contrasts sharply with the description of the connection via the Newman- 
Penrose spin coefficients [36, 44], which requires numerous differential constraint equations, 
expressing the property that the connection is torsion-free. 

The converse construction works as follows. Given /? and the connection variables {H, J, K) 
on Mz , we reconstruct u via the relation 

2 - div/3 

u = , (15) 

and we find ^,7 by solving an elliptic system for 7, 

■= + 2^ ^^"^ = ^2:^ - ^' 

and setting 

J - div7 
" = ^3dl^- 

The system (^) is M-linear and elliptic with 6-dimensional kernel, provided S/3 is not too large 
(|S/3| < (2 — div/3)/\/3 is sufficient). Prescribing the I = 1 spherical harmonic coefficients of 7 (for 
example, by requiring 7^=1 = 0) suffices to determine the solution 7 uniquely. The remaining 
connection parameter Q, together with the now known values of /3, 7 on A4, determines the 
evolution equation 

f = |^ + ^(Q + V,/?-V,7-7). (18) 
To summarise, given the field /? on a single level set A^, satisfying the size constraint 

19^1 < (2-div/3)/V3, (19) 

the map (n, v, 7, 13z) 1— > (i?, J-, K, Q, 7^=1) is invertible, assuming all fields are sufficiently smooth. 
In section |4.4| we will describe the numerical implementation of the inverse map. 
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2.4 NQS Einstein equations 

To compute the curvature of dsf^Qg, and thereby to determine the NQS form of the Einstein 
equations, we introduce the complex null vector frame {i,n,m,fh), 

n = u-^{d^-r-^-i -v{dr-r-^(5)), (20) 

m = — —{^^ — \csc{^^^p)^ 
r\/2 

and the directional derivative operators 

Vr = dr- r~^Vf3 , V, = d,- r-^V^. (21) 



Expressions for the Newman-Penrose spin coefficients [36| with respect to the frame {£,n,m,m) 
are given in terms of H, J,K,Q in Q . 

The frame components of the Einstein tensor Gab, a,b = i, n, m, fh, may be written in terms 
of the NQS metric functions and NQS connection variables. These expressions may be grouped 
into hypersurface equations (or main equations |5l[|): 

rV.H = (Ui,p.ml±I^]„, (22) 



^2 2-div/3 / 

rVrQ- = {np-uH)Q- +Q-3/3 + 2mf3 + udH - mu + 2r'^Gem, (23) 

rVrJ = -{l-dw(3)J + u-lu\Q+\'^ -ludw{Q^) -ur^Gen, (24) 

rVrK = (idiv/3 + i curl/3) K - ^3/3 + ^uQQ+ + lu{Q+f + ^ur^Gmm, (25) 

the boundary equations (or subsidiary equations) 

rV^{J/u) = v'^ rVr{J/{uv)) + \{(i].Y-f -vd].Yf3)J/u + 2u'^\K\^ 

- ^Q+v - At; + nr^Gnn, (26) 
rV^Q+ = {vrVr + J -vQ^ + Q^)Q+ - KQ+ + 2^K + 2u-^rVr{uQv) 

- (2 + icurl/3)gT; + 9 J - 2ti" VSii - 2ur'^Gnm, (27) 

and the trivial equation 

ur'^Gmrh = rVrJ -\di-vf3J -u\Q+\^ + \ud\wQ+ + Q+du + Q+du 

+ + K^() + r^{vVr - V,){u-^Vru) + u-^r^Vr{uVrv). (28) 

Observe that the hypersurface equations (^-25) have no explicit z-derivatives, and they 
each contain only one radial (r) derivative. The form of the connection variables {H, J, K, Q) 
was determined by exactly these properties. Consequently, the hypersurface equations may be 
written schematically in terms ofU = {H, Q~ , J, K) in the "characteristic-transport" form 

r^ = FiP{z,r),U{z,r)), (29) 

by treating angular derivatives such as (5U as determined by the set of values U{z,r) on the full 
(5*2. This system has the effect of transporting the fields U along the characteristic curves with 
tangent vector i which foliate the null hypersurfaces A/^. 
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Note that alternative reformulations of the hypersurface equations are possible, preserving 
the general characteristic transport structure. For reasons associated with reliably capturing 
the asymptotic behaviour of the fields, the present version of the code integrates the following 
radial equations, for the variables log(///2) (instead of H), j = jr{l — HJ) — M (instead of J) 
and rQ~^ (instead of Q~): 

rdr\og{H/2) = V^log(i7/2) + idiv/3-^^^^±^, (30) 

rdrirQ+) = Vp{rQ+)- {l-2(5p){rQ+) + V,Q+P + 2rn{rVrlogu) 

+ 2rp-[ rScurl/? + 2r^Gim, (31) 
rdrj = V isj + {j + M — ^r){dw(] — rVrlogu) 

+ lri\Q+\^ + divQ+) + lr^Gen, (32) 
2rdrK = 2Vf3K + {dw(3 + 2icuvip)K - Jdp 

+ unQ+ + lu{Q+)'^ + ur^G^„^. (33) 

Here M = 1 is a constant which fixes the bare mass of the background Schwarzschild black hole. 
Of course, the Einstein tensor components in these formulae are set to zero for the vacuum 
equations. 



It is remarkable that the boundary equations (26,27) and the trivial equation (|28| ) may be 
regarded as compatibility relations, by virtue of the conservation (contracted Bianchi) identity 
G^^'^ = |51, ^. This identity is valid for any Einstein tensor Gab, regardless of the metric. 



Substituting the hypersurface equations Gu = Gem = Gen = Gmm = into the conservation 
identity, yields equations HGmm = and a propagation system for G„„ , Gnm which has the 
unique solution G„„ = Gnm = if the boundary equations are satisfied on one hypersurface 
transverse to the outgoing null surfaces Mz- Thus in order to construct a solution of the full 
vacuum Einstein equations, it suffices to satisfy the hypersurface equations everywhere, and the 
boundary equations just on the boundary surface r = vq (for example). 



3 Numerical techniques 

In this section we describe the data representation and manipulation techniques. These consist 
mainly of techniques for handling angular fields and derivatives, and an unusual convolution 
spline used for interpolation, differentiation and high frequency filtering in the radial and time 
directions. 



3.1 Fields on S'^ 



The evolution algorithm treats the angular derivatives ^ , ^ as "lower order" , compared to the 
radial and time derivatives This attitude in a numerical computation can be justified 

only if it is possible to easily and accurately compute and manipulate angular derivatives. This 
is achieved by using spectral representations (both Fourier and spherical harmonic) for fields 
on the 2-sphere. This approach is widely used in geophysical and meteorology applications 
P^ , |4^ , |l5| , [5^ , 61 1 and is known to have significant advantages compared to finite difference 
approaches ||T^, based on either angular coordinate grids or overlapping stereographic projection 
Nevertheless, spectral methods have rarely been used in numerical general 




charts p5 

relativity (however, see |45, 13 1) and they have not been used previously for solving the full 
Einstein equations. 
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The basic manipulations required of S'^ fields are: 

• computing non-linear algebraic terms such as f/(2 — div/3), ulQ"*"!^ etc; 

• computing angular derivatives operators such as div/3, (5Q~^ etc; 

• inverting the linear elliptic operator £^ (which appears in equation ([l6|)); and 

• projecting aliased or noisy field value data onto certain subspaces of spin- weighted spherical 
harmonics. 

To carry out these manipulations, three separate representations are used for fields on S'^: 

• field values r]jk = ri[{}j,ipk) at the polar coordinate grid points 

(^?,-,(^fc) = ((j-i)AT?,(A:-l)A(^), (34) 

where At? = 2tt/N, Aip = 2Tr/N with 1 < j < N/2, 1 < k < N and (in our implementa- 
tions) N = 16, 32 or 64; 

• Fourier coefficients r]mn arising from FFT transforms in the •& or (p directions, of the field 
values rjjk] 

• spin-weighted spherical harmonic coefficients rf '"^, \m\ < I, I = s, . . . , L, L = N/2 — 1 (for 
spin s = 0, 1 and 2). 

The field value representation is used when computing non-linear algebraic terms such as 
ulQ"*"!^. The Fourier representation is used for computing -i? and ip angular derivatives, which are 
needed in the formulas for S, div, for example. The spherical harmonic representation is used 
in solving the elliptic system (p!6|), to spectrally limit the field values by projection to spherical 
harmonic data, and to summarise the computation results (which are stored using spherical 
harmonic coefficients). 

For fields which do not alias on the (i9, ip)-gnd, the three representations are completely 
equivalent in the sense that conversion between them is essentially exact, depending on the 
machine precision and on algebraic details of the specific FFT algorithm used. The requirement 
that a field does not alias is satisfied when it can be represented by a finite expansion in spin- 
weighted spherical harmonics with angular momentum I < L = N/2 — 1. Our implementations 
use spectral cutoffs L = 15 (for N = 32) and L = 31 (for N = 64). 

Transformation to the spherical harmonic representation involves a projection, because both 
the field value and Fourier representations have approximately twice as many degrees of freedom. 
For example, a (real) spin field with I < L = N/2 — 1 has (L -|- 1)^ = N'^/4 spherical harmonic 
coefficients, whereas it has N'^/2 values on the (??, (/3)-grid. The space of non-aliasing spherical 
harmonics is a linear subspace of the space of functions represented by either Fourier coefficients 
or field values. 

For example, when the (■(?, ip)-gnd field values of the product of two fields is calculated, the 
result, which will contain components up to / < 2L, is aliased onto the grid in such a way that 
its field values no longer lie in the appropriate spin weighted spherical harmonic subspace. To 
clean up after such non-linear effects, we project the result back onto the correct subspace, as 
described in section |3.5| . 

To minimise the possibility of unstable feedback of quadratic aliasing errors, we may invoke 
the Orszag 2/3 rule |jl^, ^ at various points within the code. The effective spectral resolution 
of the code is then /max ~ 2L/3 (/max = 10 for = 32 and /max = 20 for N = 64). 
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3.2 Spherical harmonics 

We first summarise the more important properties of 9 ("edth") and spin weighted spherical 
harmonics. The edth formalism provides a unified geometric approach to the treatment of 
angular derivatives on 5^ and vector and higher-rank tensor harmonics. Detailed descriptions 
of the properties of spin-weighted fields and spherical harmonics may be found in Penrose and 
Rindler |44| or [23, Here we describe only the basic formulae. 

We use a real-valued basis Yim, Z = 0, 1, 2, . . ., m = — /, . . . , Z for the space of spin spherical 
harmonic functions, defined by 



where 



1 m = 

cos mip m > 
■v/2sin|m|99 m < 0, 



(35) 
(36) 



and the Pim{'&) = Pi\m\{'&) ^'^^ related to the associated Legendre functions Pim by 



l\m.\ 

Plm{^) 
PlmicOSt)) 



(-i)™^/2rTi 



'{l-m)\ 
{I + m)] 



2H\ 



■ sin™ ^ 



[X 



Pim (cost?) , 



1)' 



The spin s spherical harmonics are then defined explicitly by 



Y, 



Im 



(-ir 



2'{l-s)\ 
2'{l-s)\ 



1/2 
1/2 



Q'Y, 



Im ■ 



s > 0, 
-s < 0, 



(37) 
(38) 

(39) 
(40) 



where necessarily / > Note that the differential operator S is spin-raising, sending spin s 
into spin (s + 1) fields, and S is spin-lowering. 



[i(Z + . + l)(/-.)]^/^l^^+\ 



S^/m — ~ L2 V T •'A' ~ "T -l;J ^Im ' 

for all s G Z, and Yj^ and Yj^ are related by complex conjugation. 

Since Al/^ = + 1)1/^) the fundamental commutation relation 

[3, S]r/ = (8S — S9)r/ = sr], 
for any spin s field rj, may be used to show that 

^YC^ = {s^-l{l + l))YC^, 
where A = S9 -|- 8S. With these conventions we have the orthogonality relations 



(41) 
(42) 

(43) 

(44) 

(45) 



47r 



T- / Yi'm' sin ^d'dd^ = 8iv5„ 



52 
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which show that the YJ^ form a basis (over C) of the Hilbert space of square-integrable spin s 
fields on 5^, which is orthonormal in the natural Hermitian inner product 

{(P,^) = ^i (46) 
J s 

Prom (|3|-||) it is evident that the spin harmonics Yim are trigonometric polynomials in "d 
and (p [34|. Using expression (^) for S, we also see that the are trigonometric polynomials. 
The highest wave number Fourier modes which occur in the set of basis functions {Yf^ '■ s < I < 
L, \m\ < 1} are cos(L'i?), sin(L??), cos(L(^), and sin(L99). Therefore, on a uniform {'&,ip)-giid 
of size N/2 x N, one can represent all of the spin weighted spherical harmonic functions up to 
L = N/2 - 1. 

3.3 Even/Odd decomposition 



Because 8 is surjective onto the space of smooth spin s fields for s > 1 [44|, the decomposition of 
spin fields or functions into real and imaginary parts may be propagated to higher spin. The 
resulting decomposition into even and odd components plays an important role in the analysis 



of the linearised Einstein equations about the Schwarzschild spacetime [47|. Because the NQS 
geometry distinguishes the Schwarzschild metric and is also based on spherical harmonics, it is 
ideally suited to comparing nonlinear evolution to the comparatively well-understood black hole 
linearised Einstein equations ||l^, It is thus not surprising that the even-odd decomposition 
proves to be very important in analysing the results of the NQS evolution. 

We say that the spin s field r] is even if rj = Q^f for some real- valued function /, and r/ is 
odd if f] = iS^g for some real-valued function g. (If s < then we interpret as (— S)'*'). This 
matches the usage in [47| — note that for axially symmetric fields the terms polar (for even) 
and axial (for odd) are sometimes used The surjectivity of S onto spin s > 1 ensures that 
every spin s field may be uniquely decomposed into a sum of even and odd parts. For s = 1 this 
decomposition corresponds exactly to the classical Hodge-Helmholtz decomposition of a vector 
field into the sum of a gradient and a dual gradient (or curl) — see Table H for a summary of 
the various nomenclatures. 



Even: 


polar 


irrotational 


9/ 


grad/ 


{Vlf)vi + iV2f)v2 


Odd: 


axial 


divergence-free 


idg 


curl(/ 


(V25) vi - {Vig)v2 



Table 1: Equivalent terminologies for vector fields on S'^ 

The even/odd decomposition has a natural interpretation in terms of the spectral decompo- 
sition 



l=s m=—l 



of a spin s field r/, because we are using a basis of real- valued YJ^. Namely, r] is even if the 
spectral coefficients rf"^ are real, and odd if the coefficients are pure imaginary. We sometimes 
use Even(r/) and Odd(77) to represent the respective projections, so r/ = Even(77) -|- Odd(77) with 

Even(r?) = J^Re(r/'-)y4 (48) 
Odd(ry) = i^Im(V™)y4. (49) 
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Observe that QY^^ = for s > and thus S acting on spin s fields with s > has kernel 
having (complex) dimension 2s+l. Likewise the formal adjoint —8 acting on spin s < fields has 
(2|s| + l)-dimensional kernel. In particular, S acting on spin 1 fields has kernel consisting of the 
C-linear space spanned by the three I = 1 spin 1 spherical harmonics — the corresponding 
real vector fields are the dual gradients of functions linear in M^, which are just the infinitesimal 
rotations, and the gradients, which are the conformal dilation vector fields. 

The correspondence between vector fields on S'^ and spin 1 fields generalises to spin 2 fields, 
which correspond to symmetric traceless 2-tensors on S'^. If A is a symmetric traceless 2-tensor 
then with respect to the standard polar coordinate derived orthonormal frame 



we have the correspondence 



ei = d^, 62 = csc&d^, (50) 



i(Aii-A22-2iAi2). (51) 



This correspondence extends to higher integer spins with higher rank symmetric traceless tensors 
on S'^. The cases s = 0,1,2 of most importance in our work correspond to the usual scalar, 
vector and tensor harmonics. 

3.4 Fourier representation 

The fast Fourier transform (FFT) is used to transform between Fourier coefficients and field 
values on the uniform ("i?, ip)-gnd. Fourier convergence problems arising from discontinuities in 
coordinate derivatives and vector and tensor components at the poles, are sidestepped by an 
observation relating fields on 5^ to fields on the torus = 5^ x 5^ [^]. The torus method 
enables coordinate and covariant derivatives for all types of field to be computed using Fourier 
methods, so is particularly well-suited to handling the derivative operator S (^). This approach 
to handling component discontinuities at the poles is simpler than the techniques reviewed in 



1 59 1 for manipulating vector fields, and readily extends to any rank s > 0. 

For integer s the real and imaginary parts of a spin s field on S'^ may be identified with the 
two independent frame components of a completely symmetric trace- free tensor of rank |s| on 5^ 
|44|. Since the frame ei,e2 (^) is not continuous at the poles, the tensor components will not 



be continuous at the poles, so are not obviously suited to Fourier expansion in the i? direction. 

However, along any smooth curve crossing through a pole, both basis vectors ei and 62 
reverse direction at the pole. Thus, for any smooth tensor field T = T^^'^^^'Cj^ • • • ej^, by 
continuity of T the component functions T^^---^" will change by a factor (— 1)'^ across the poles. 
Consequently, if we extend the domain of definition of T-^i - -?' to € [— vr, vr] by 

TJi-J'-(-i9,(^) = (-l)^r^'i-J'-(^,^ + 7r), fort?G[0,^] (52) 

(using the 2tt periodicity in ip), then the resulting extension is 27r-periodic and continuous in -d. 
This argument extends to higher (covariant) derivatives of T, showing that the extension is in 
fact smooth and periodic in ■!?. Derivatives of T-^i----?^ with respect to can then be calculated just 
as for if derivatives, provided that the direction of increasing is properly taken into account. 
In effect, the extension just described defines a (smooth) field T-'i on the torus = 
X S^. This may be understood geometrically by noting that the map 
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is in fact smooth. This follows by noting that because is a radial coordinate near the north 
pole "i? = 0, the differential structure near the pole is represented by the rectangular coordinates 
(^,77) = (i? cos "i? sin and the map 93) 1— > {(,,rj) is manifestly C°° for i? near 0. Conse- 
quently any smooth tensor T on S^, when expressed in a cotangent basis, pulls back to a smooth 
tensor on (ie. T*(r) G C°°(T^)), and thus admits a well-behaved Fourier representation on 
T^. The converse is of course false: a smooth tensor on does not necessarily arise from a 
smooth tensor on S*^, even if it satisfies the parity condition ( |52D satisfied by pull-back tensors. 

The coordinate derivative form (^) of S (and similarly any other covariant angular derivative 
operator) can be easily evaluated by transforming to the Fourier representation of the field, 
multiplying the Fourier coefficients by the appropriate wavenumber factors, transforming back 
to obtain the field value representation of the "d and (p derivatives, then finally including the 
cscd factors. Thus, computing a derivative operator such as 9 is an O(A^^logA^) operation. 
Since csc(iA??) ~ 10 for N = 32, there is no significant loss of accuracy in calculations near the 
poles. 

It should be emphasised that because the are trigonometric polynomials in {'d,(p), the 
FFT computation of their numerical derivatives is algebraically exact. For example, the Lapla- 
cian relation (|45| ) is numerically verified to 1 part in 10^^ |39| ]. 



3.5 The spherical harmonic representation 

The field value and Fourier representations suffice for most numerical calculations, such as 



computing the nonlinear terms in (22-25). However, at some steps it is essential to use the 



representation by spherical harmonic coefficients (|47]): 

1. to solve the elliptic system (1T6|) (see section |4.4| ); 

2. to implement spectral projections after each time step, with the aim of suppressing aliasing 
effects and rounding errors at the poles ; 

3. to store results for later analysis and display, since the spherical harmonic representation 
is more compact and the coefficients may be readily interpreted physically, by comparison 



with well-studied solutions of the linearised black hole Einstein equations [47, 65, 19 1. 



We next describe the methods used to transform fields from the field value and Fourier 
representations to the spherical harmonic representation. Transformations for fields of spin 0, 1 
and 2 are required by the code; the spin case which we describe here to illustrate the technique 
is slightly less complicated, since we may assume the field / is real-valued. 

There are (L + 1)^ basis functions in the set {Yi^ : < \m\ < I < L}. However, to 
represent these functions as trigonometric polynomials on a regular S'^ grid we require a grid 
of size (L + 1) X 2(L + 1), and thus 2(L + 1)^ real coefficients. The spin functions of angular 
momentum at most L therefore form a subspace of real dimension (L + 1)^ in the space of Fourier 
series representable on the grid, which has real dimension 2(L + 1)^. We use a projection onto 
the spherical harmonic subspace which is orthogonal with respect to the natural inner product 
in the Fourier space, 



1 



2iT r2iT 



(/l,/2)2 = ^/ / fl{'&,V)f2{'&,V)dM^ 



JO 

N N 



]^EE/i(^-'^.-)/2(^.,^,), (54) 

2=1 j = l 



12 



where ipj) : i, j = I, . . . , N} are grid points (cf. (||)) and N = 2{L + 1). 

To make use of ( |5^ we use (52) to extend functions defined on to functions defined on the 
torus = S^xS^. In particular, given any set of values {fij G M : i = 1, . . . , N/2, j = 1, . . . , N} 
on the S"^ grid, we use (52) to construct grid values on T^. There is then a unique interpolating 
trigonometric polynomial / such that /(i?j,(/9j) = fij. We project / to the / < L spherical 
harmonic subspace as follows. 

The are not orthonormal with respect to (p4[), but instead have Fourier inner product 



Glml'm' — (XlnnYl'm')Y 

= (PlM,Pl'mmYSmm'- (55) 

Here, the index pair Im (and I'm') is a combined index which takes {L + 1)^ values, and (|55| ) is 
the matrix for the induced Fourier metric on the spin subspace. The summation convention 
will be employed for raised and lowered repeated indices. 

For fixed m, the inner product (|55|) of the P functions forms a matrix. 

These matrices are defined only for \m\ < 1,1' < L, so are square and of size (L + 1 — |m|) x 
(L + 1 — \m\). Denoting the inverse matrix by we have 

Glm I'm' — A{m)ll' '-'mm' ) 

(56) 

and the components of the inverse metric are given by 

Qlm I'm' ^ ^11^^ ^mm' (^^^ ^) ^ (57) 

The dual basis vectors for the spin subspace are 

y'™ = G^"''™'y«/„/, (58) 

and satisfy (l/m, )f = ^m ■ The orthogonal projection of / onto the subspace is given 
by 

proj(/) = (/, Y'"')pYi^ = f^Yim, 

where 

fm ^ Y'm^p (59) 

are the spherical harmonic coefficients of the function /. 

To calculate the inner product (|5^), first note that using (|35|), (|57|) and (^), the dual basis 
vectors can be written as 

y'™ = p'"'(i9)F'"((/;) (60) 
(in analogy with (|35|)), where we have set 



P'™(^) = Af;^)Pi'mm , F"^{^) = FM. (61) 

By Fourier analysis of / in the direction we can write / = f^{'d)Fi.{(p). In particular, by 
(y3-FFT of {/jj} we get the numbers f^{'&i)- The spectral coefficients can then be evaluated 
using ( p4[ ) and ( |60| ) as 

fim ^ (/fc(^)Ffc(v9),p'"^(^)F-(<^))F 

= (/-(^), p""(^))f 

= (62) 
1=1 
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The converse process of reconstructing the function values fij = fl'&ijipj) from the spherical 
harmonic coefficients Z'™" follows from 

= f^PiMFM- 

First we construct the quantities 

L 

^{^^) = Yl ' (no sum on m), (63) 

l=\m\ 

and then we use inverse FFTs in the (p direction to reconstruct fij via 

L 



in=—L 



Both constructions, of J''" from fij and conversely, are 0{L^) operations, due to the matrix 
multiplications in (^), (|63|). Routines for transforming between grid values and spherical har- 



monic coefficients have been implemented for maximum angular momentum L = 7, 15 and 31. 
The grid values P {"di) which appear in the sum (^) were pre-computed in multiple precision 
using REDUCE |Q. The functions P defined by (|6T|), were constructed symbolically using 
exact inversion of the matrices This symbolic approach was feasible because the metric 

Gim I'm' factorized as the tensor product (|55|), thus allowing exact inversion of G using matrices 
of size at most (L + 1) x (L + 1) rather than {L + 1)^ x (L + l)^. 

The analysis of spin 1 and spin 2 grid functions into spherical harmonic coefficients is similar, 
but complicated by the fact that the induced metric on the subspace factorizes as a tensor 
product only in a complex (mixed parity) basis. Separating the even and odd parity coefficients 
therefore requires some extra book keeping. 

Techniques for handling spherical harmonic spectral representations have been described by 



many authors |34, 42, 17, Our method differs from the Muchenhauer and Daly projection 
(see [58 1) in the choice of inner product (U) used to define the orthogonal subspace. More general 
spectral transform methods (eg. |T^) use other choices of weightings and node points to define 
the projection, and do not have such a simple underlying inner product. All these methods are 
also O(L^). Jakob [^] gives an O(L^logL) spectral projection, which however bypasses the 
construction of the spherical harmonic coefficients. Since we need the spectral coefficients, and 
because we work with a relatively small value of L, the Jakob projection would not provide any 
improvement. 



The torus method described here and in |3£] has the advantage that it applies also to higher 



rank tensors, in particular vectors and 2-tensors. Representations in terms of spin-weighted 



fields are more efficient for vectors (spin s = 1) than 3- vector representations |59, 60 1, and the 



operator 9 gives a transparent derivation of all invariant derivative combinations |59]. 



3.6 Convolution splines 

At various stages it is necessary to interpolate and differentiate grid-based fields. For example, 
the radial integration of the hypersurface equations by the 8th order Runge-Kutta method 
requires values of the source field /? at 10 intermediate points; the dynamic regridding of the 
radial grid requires interpolation to determine the field values of P at the new grid points; 
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and derivatives such as d^/dr and dQ^ jdz must be computed from field values on grids. A 
convolution spline algorithm described in provides a convenient technique. 

The method has the effect of fitting a spline curve to sample data, and is implemented by a 
convolution of the form ||3^ 

f{x) = Y,f{k)K{x-k), (64) 

fcez 

where the /(A;) are the raw data (samples) and 4>n{x) is a (7""^ sampling kernel. 

The sampling kernel 0„ is constructed as a certain sum of central B-splines M„ of order n, 

n-l 

Mx) = J2 «S"^^n(x - t + i), (65) 

i=l 

where the coefficients a[^\ i = 1, . . . ,n — 1 are chosen so that the convolution ( |64[ ) acts as the 
identity on polynomials f{x) of degree n — 1 (n even) or n — 2 (n odd). Recall that the central B- 
spline Mn{x) is a C"~^ piecewise polynomial of degree n — 1 normalised by Ylik&'L ^n{x — k) = l, 



with support on |x| < n/2 |52]. The support of the kernel (j)n{x) is therefore |x| < n — 1. 



Algorithms for computing the a^"^ and tabulations for n < 11 are given in ||38|]. Coefficients 
for the kernel (j)g used in the code are given in Table ^, and (pg is plotted in Figure [l|. 



Table 2: Convolution coefficients a\ 



i: 


1,8 


2,7 


3,6 


4,5 


(9) 


67 
2520 


nil 

5040 


421 
560 


1333 
1260 



Sampling kernels: 0g(x) and sinc(x) 




-5 5 

X 

Figure 1: Comparison of the spline kernel (j)g (solid curve) and the sine method |[l^, |5^ kernel 
sinc(a;) = sin(7rx)/-7rx (dotted curve). The kernel for a sampling method is the response to the 
delta-like discrete data (solid dots). The advantages of the convolution spline method are that 
the kernels have finite support ([—8,8] in this case) and the data is automatically filtered. The 
convolution (^) using kernel (pg exactly reproduces polynomials of degree 7 or less. 

The expressions (|6^), (^5|) may be rearranged into a form which is more efficient for numerical 
calculations, 

f{x) = Y,fkMn{x-^-k), (66) 

fcez 
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where the modified sample values are given by 

n-l 



i=l 



The advantage of (|66[) is that the fk can be computed once and then reused to evaluate /(x) at 
many different points x, using the explicitly known values of the B-spline M„(x). The same fk 
values may also be used to compute derivatives of the spline function /, 

P\x) = Y,fkM^^\x-k-^), (67) 

fcez 

where again the derivatives Mn\x) are known functions. These techniques are routinely used 
to supply intermediate values and derivatives of fields in the radial and time directions. 

Non-uniform distributions of sample points are handled by a mapping between the indepen- 
dent variable and the sample number variable (x in the above). Numerical derivatives are then 



calculated using the chain rule. For example, the radial grid described in Section 4.2 is non- 
uniform, specified by some known relation of the general form r = r[n) (where n is now being 
used to denote the sample number variable, with radial grid points being at n = 0, . . . ,noo)- 
The operator ^ is then implemented as (^) with a formula for the first factor being 

known explicitly. 

Similarly, one can transform to an independent variable s = SQ + hx which has grid spacing /i, 
to examine the behaviour of the approximation (|6^) as /i — > 0. Let g{s) := f{{s — so)/h) = /(x), 
so g^-'\s) = f^^\x). It can be shown ||3^ that using the (j)g kernel, the Taylor series truncation 
errors for ( |6^ at a grid point s are 

\-gU)(^s)-g^^\s)\ = cjhV'Hs)\+0{h^), j = 0,1,2, (68) 

where cq = jlffgQ, ci = C2 = ^^q- This reflects the fact that convolution with (f>g is 

exact on polynomials of degree 7. 

The predicted convergence of the (jyg spline convolution is clearly evident in Figure |2|. Here 
the function v{x) = sin lOx has been approximated at varying grid resolutions corresponding 
to = 2^ grid points over the interval [—1,1]. 

Convolution splines do not generally preserve sample values, except for samples from poly- 
nomials of degree less than or equal to the degree of reproduction. This results in some damping 
of high frequency components of the data, which we expect helps to suppress numerical noise 
and algorithmic instabilities. 

Within the context of spectral methods for PDEs, the direct filtering of Fourier coefficients 
of a numerical solution is common practice and has been extensively studied (cf. |l8| , §8.3] and 
references). On the other hand, explicit use of a digital filter [^] in conjunction with finite 
difference methods is comparatively rare. Nevertheless, from an algorithmic point of view, this 
is the effect of using a convolution spline. 

The filtering inherent in the method can be examined via the response function 

^.„(6') = V (l),,{k) cos ke, 0<e <Tr, 



which is equal to the factor by which the Fourier mode e'^^ is amplified by the approximation 



(64) at a grid point x G Z. The value 9 = ir corresponds to the Nyquist frequency for the grid. 
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V = e''sin(10x), N=16 



V error for N = 2'', p=4 7 




-1 -0.5 0.5 1 1-1 -0.5 0.5 1 



X X 

(a) (b) 

Figure 2: Convergence of the 4)g spline convolution: (a) samples of v{x) = e^sinlOx at A'^ = 
16 points over [—1,1], and the corresponding convolution spline; (b) logarithmic plots of the 
absolute error \v{x) — v{x)\ for grid resolutions of = 2^ with p = 4, . . . , 7, showing a reduction 
in the error by a factor of 2^ on each doubling of the resolution. 

Filter response functions 




0.2 0,4 0.6 0.8 

frequency, S/tt 



Figure 3: Filter response functions: (a) raised cosine filter (j{0) = 2(l + cos^); (b) Lanczos filter 
a{9) = O^^smO; (c) sharpened raised cosine [^]; (d) ^<i{9). 

Figure ^ shows the response function <&9(0), compared to the well-known Lanczos and raised 
cosine (artificial viscosity) filters p9| . The filtering characteristics of convolution splines and 
their derivatives are described in [p^ ]. 

The limitations of convolution splines are illustrated in Figure ^ which shows Gibbs-like 
effects associated with approximation of step-like data. The figures also give an indication of 
the number of grid points needed to resolve a sharp transition in field values. 

In order that the convolution spline smoothing should introduce only negligible errors, the 
grid resolution should be chosen sufficiently fine that typical field variations take place over 
enough grid points that the expected frequency 9 of the field data lies well within the part of the 
Nyquist frequency interval where ^n{9) ~ 1- Of course this requires some prior knowledge of 
the length scale of the field, and cannot be applied where shocks (or arbitrarily rapid variations) 
occur in the data. In such cases the spherical harmonic representation would become equally 
unsuitable. 

The choice of high order convolution splines (/i^ rather than say h^) was motivated by the 
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u = tanh(IOx), N = 16 



N = 64 




Figure 4: Convolution spline approximation of step-like data. Here the function u = tanh(lOx) 
is represented by N samples over the interval [—3,3]. For N = 16, the transition from n ~ — 1 
to 1 takes just one grid interval and a Gibbs-like phenomenon is evident. The convolution 
splines are constructed using the 4>9{x) kernel. The logarithmic plot shows the absolute error, 
\u(x) — u{x)\. 



need to reduce storage requirements. Low order spline convolutions have a markedly reduced 



usable proportion of the Nyquist interval [38|, so to avoid over-smoothing and to achieve compa- 
rable accuracy with a low order method would require significantly higher grid point densities. 
Both storage costs and the cost of the radial integration increase linearly with the number of 
radial grid points. The (pg kernel was chosen also because its accuracy matches that of the RK8 
method used for the radial integration. 

To use convolution splines near endpoints of a data set, one can extend the data set, using a 
suitable mapping between the independent variable and the sample number variable |^^. The 
mapping is chosen so that when expressed in terms of the sample number variable n, < n < Uoo, 
the fields admit expansions in powers of near n = and (noo — n)'^ near n = n^. The 
sampling kernel convolution can then be applied to the even extension of the fields through 
n = or n = noo- This technique is particularly important in extracting radiation data near 
null infinity (scri, X"*"), where the radial grid is chosen so that noo — n = 0(r^^/^), as described 
in Section ]A:.2[ 



4 Solution algorithm 



The hypersurface equations (p^-p^) suggest the following process for evolving the metric in the 
exterior region with interior boundary the cylinder r = tq: 
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1. Choose boundary data {H, Q , J, K) on the cyhnder r = tq, consistent with the boundary 
equations (2^,27[); 



2. Assume P is given on a null hyper surf ace J\fz', 

3. Solve the A4 hypersurface equations rdrlJ = F{I3,U), by integrating along the radial 
curves (z, = const, with initial conditions at r = rg determined in step 1; 

4. reconstruct the metric functions u, v, 7 from H, J, K and [i using the converse construction 




5. reconstruct 8(3/ dz from Q and the now known values of /?, 7 on Mz^ using (p^); 

6. use 813/ dz from step 5 to evolve 13 to the "next" null hypersurface Mz+Az and repeat from 
step 3. 

In the following we will show how this heuristic algorithm is implemented numerically, using the 
techniques and data representations of the previous section. 



4.1 Geometry and inner boundary conditions 

The code models gravitational waves propagating on a black hole spacetime, with metric approx- 
imating that of the Schwarzschild solution in the Kruskal-Szekeres coordinates ||30t| . Introducing 
the double-null coordinates 

z = t-r*, y = t + r\ r* =r + 2Mlog - 1) , 

the Schwarzschild metric becomes 

dsl^h^ = -(1 - 2M/r) dydz + r'^dQ^ , 

where dO,"^ = dd"^ + sin^ -d dip"^ . The coordinate singularities at the past and future horizons 
t = ±00 are removed by defining y = e^/^^, z = e~^/^^, giving the metric 

dsl^^^ = ^e-^/2M dydz + dQ^ , 
where r = r{y, z) is defined implicitly by (see Figure ^) 

.-'■^"(i-O-P. (TO) 

Because radial light rays are straight lines at 45° and r = is singular, the surfaces 5 = and 
y = (ie. r = 2M) form the past and future event horizons, and these are smooth hypersurfaces 
with bounded curvature. The approximate Minkowski structure of Schwarzschild spacetime is 
better illustrated by the radial null geodesies in (r, t) coordinates, see Figure |6|. Note however 
that the (r, coordinates are singular along the event horizons r = 2M. 

Initial conditions for (3 are imposed on {z = 0,r > 2M} (with M = 1 usually), by specifying 
the spherical harmonic coefficient functions (3im{r). Since (3{z = 0) is unconstrained, these 
coefficient functions may be freely chosen, subject only to the size condition (p!9[). 

For simplicity the inner boundary conditions are set at tq = 2M = 2 to agree with the 
Schwarzschild past horizon: Hq = 2, Qq = Jq = Kq = 0. Since we choose /3(0, r) = for 2 < r < 
5, by causality the solution should agree with the Schwarzschild metric in a neighbourhood of 
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Kruskal— Szekeres coordinates 




Figure 5: Schwarzschild spacetime in Kruskal-Szekeres coordinates. Radial light rays are straight 
lines at 45°. 



r = 2 for all time z > 0, producing a "white hole" past horizon in the spacetime. This choice of 



inner boundary condition has the considerable advantage that the boundary equations (|26D , (|27|) 
are automatically satisfied, and it is not necessary for this class of simulations to separately 
ensure that the boundary data are numerically compatible with the boundary equations. 

The resulting spacetimes have the geometry of an isolated black hole with future event 
horizon at z = oo, r > 2M and Schwarzschild-like white hole boundary along r = 2M,0 < 
z < oo. Adding shear (3 at the initial hypersurface z = results in spacetimes modelling the 
interaction of gravitational radiation with a single black hole. 

Although the fixed past horizon boundary conditions used in the present code allow many 
interesting issues to be addressed, it would be desirable to implement more general inner bound- 
ary conditions. Such conditions specify H, J, K,Q~^ at an inner surface (r = 1, for example), 
subject to the dynamical {d/dz) constraints on the evolution of J/u and determined by 
the boundary equations (p6|),(p7|) Q. The free data on the inner boundary consist of Uo,Kq, 
where uq represents a certain coordinate gauge freedom, whilst Kq describes the gravitational 
radiation injected into the system through the inner boundary. Various exact solutions with 
such boundary conditions are described in Q (Robinson- Trautman |49], boosted Schwarzschild, 
twisted Minkowski space), and would provide useful accuracy checks on the numerical methods. 

However, implementing general inner boundary conditions raises numerical and geometric 
difficulties — the boundary data must be "consistent" with the RK4 solution evolution in order 
to maintain optimal accuracy (see for an analysis of similar but simpler situations), and 
constraining the radiation data Kq such that the spacetime is still Schwarzschild near the past 
horizon is a difficult geometric problem. An arbitrary choice of Kq (even Kq = if uq 7^ 1) will 
inject some additional energy into the spacetime. 



4.2 Dynamic radial grid 

There are two geometric features which the code should model accurately: future null infinity 
("scri" or 2~^, where r — > 00, z finite), and the future horizon r ~ 2M, z ^ 00. 

The field near null infinity X+nA4 = (r = 00, z) determines the outgoing gravitational waves 
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Schwarzschild coordinates Finkelstein diagram for (r,z) coords 




(a) (b) 

Figure 6: Exterior region r > 2M of Schwarzschild spacetime, with past and future radial 
null geodesies, (a) In Schwarzschild coordinates {r,t); (b) in retarded Eddington-Finkelstein 
coordinates (r, z). 



as seen by a distant observer and is consequently very important for applications to gravitational 
wave astronomy. Experience with 3+1 codes shows that it is not possible (as yet) to provide 
boundary conditions on an outer timelike boundary at a finite radius which do not either inject 
radiation or reflect radiation back into the grid. This deficiency has the effect of severely limiting 
the overall time duration of most 3 + 1 simulations. We avoid all such reflection problems by 
using a radial grid coordinate n, which compactifies r = oo and leads to accurate modelling of 
gravitational radiation. 

The {z, r) coordinates become singular near the future horizon z — oo in the Schwarzschild 
spacetimes (see Figures ^, |6|). In our case this picture is not exact, since the spacetime geometry 
is only approximately Schwarzschild, and thus the future horizon (for example) will not be 
located exactly at z = oo. However, the NQS parameterisation must still become singular 
eventually, as the outgoing null hypersurfaces A4 approach the future event horizon. 

The effect of the nearly singular coordinates is that at late times, the in-falling gravitational 
components near the event horizon will be compressed into a region of small r-variation, and 
this compression will accelerate in time z, whilst retaining field structures from early times. 
Consequently no r-grid which is constant in time is able to accurately represent the in-falling 
radiation at late times. We have observed that numerical problems with a fixed radial grid arise 
as early as z = 10. 

To overcome these problems, a dynamic and variable radial grid is used, based on dou- 
ble null coordinates {z,y). The time steps in the evolution direction are regular, with Az = 
0.1,0.05,0.025 being typical. The grid in the radial direction is chosen to satisfy the criteria 
that it compactify null infinity and concentrate grid points in the region of greatest variation 
in the seed field /3. Because the field features propagate along the inward and outward null 
characteristics, which correspond respectively to the curves y = const, (approximately!) and 
z = const, (exactly), the numerical grid is taken to be rectangular in the {z,y) coordinates, 
with radial grid point positions being determined by an initial distribution of grid points on the 
surface z = 0. 
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Introducing the radial grid coordinate n with range < n < rioo (with typical values of rioo 
being 128, 256, 512 and grid points at integer n), we specify an initial grid point distribution 
r{z = 0,n) = f{n), where / is some monotone increasing function such that /(woo) = oo. The 
radial grid points on the initial {z = 0, z = 1) surface have y ordinates given by (^ 



y = (/(n)/2M - 1) exp(/(n)/2M) = 0(/(n)/2M), (71) 

where (p{x) := (x — l)e^ is monotone and invertible for x > 0. Since y,z,r are related by 
y = e^/^^(/>(r/2M) and the grid points are required to inflow along the curves of constant y, we 
can determine the dynamic radial grid point distribution r = r{z,n) in terms of the initial grid 
distribution function f{n) by 

r(z, n) = 2M 0-^(exp(-z/4M)(^(/(n)/2M)) , (72) 

where the inverse function : [—1, oo) [0, oo) is evaluated numerically. With this definition, 
the surfaces n = const, correspond to in- falling null hypersurfaces in the reference Schwarzschild 
metric. In order to express the hypersurface equations in terms of n rather than r, we need to 
compute dr/dn — this expression follows immediately from ([7^): 



^ _ ^-z/iM fjn) _ ( f{n)-r{z,n) \ 

dn~^ r(z,n) 2M )dn ^"^> 

It remains to choose the initial grid distribution f{n). The condition riao — n = 0(r^^/^) 
is achieved by setting f{n) = fi{v)/{l — u)'^ , where v = n/rioo and /i : [0,1] ^ M is any 
suitable smooth monotone bounded function. In the code, /i is a quadratic polynomial, with 
coefficients chosen to concentrate grid points across the support of the chosen initial data I3{z = 
0). Figure |^ shows sample curves r{z,n) = const., illustrating the in-falling nature of the 
(z, n) grid coordinates. This heuristic prescription for distributing the grid points works well in 
practice — Figure ^ shows the shear over the {z, n) plane for run_160, and clearly demonstrates 
the in-falling structure of this solution. The simulation eventually terminates at z = 55 because 
of some geometric effect associated with breakdown of the NQS gauge condition near the future 
event horizon. 

4.3 Hypersurface equations 

The hypersurface equations are solved by treating them as a large system of ordinary differential 
equations, with the radial grid coordinate n playing the role of independent variable, and the 
dependent variables being the values taken by the fields {H, J, Q, K) at the N"^ /2 points of the 
(??, (/?)-grid. 



The form (30-3^ of the hypersurface equations, for the variables logH, rQ~^, j and K, 
proves to be better behaved near r = oo, since each of these variables has a finite (usually 
non-zero) limit. Integration of these radial ODEs is possible up to and including the final point 
n = noo, with results whose numerical effectiveness may be seen by inspecting the field values 
in a neighbourhood of null infinity [^]. Tests described in the following section, in particular 
the consistency of the constraint equations and the accuracy of the Trautman-Bondi mass decay 



formula (Figure 17) also confirm that asymptotic behaviour has been reliably calculated. 

Note that unlike methods based on Bondi-Sachs or Newman-Unti coordinates p8| , |37[ ], in- 
tegration along the r-coordinate lines does not correspond to integrating along the radial null 
geodesies (the characteristics of the Einstein equations), since in general the NQS shear f3 is 
non-zero, and the null direction is £ = d/dr — r~^/3. 
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r=const. curves in the (z,n) plane 




n 



Figure 7: The map between radial grid point number n and radius r is dynamic, chosen so that 
grid points approximately follow inward null geodesies. At late null-time z, grid points cluster 
near the black hole horizon at r = 2. 



The radial integration with respect to the coordinate n is performed using an 8th order 
Runge-Kutta scheme [46|, with RK step size An = 1. The RK8 method requires 13 derivative 
evaluations per RK step, of which 10 are at intermediate points not on the radial grid. Values of 
the field f3 and its angular derivatives at these intermediate points are provided by convolution 
splines generated using the kernel (pg with samples at integer n. 



4.4 Reconstructing the metric 

Step 4 of the solution algorithm requires us to reconstruct the metric functions {u, v, 7) from 
the solution {H, J, K, Q) of the hypersurface system (p^-[25|) with seed (3 and boundary data 
(i/, J, (5)|r-=ro • Note that the connection variables (1^1^) are determined by the values of 
the metric functions {(3,^,u,v) and d(3/dz on the hypersurface N-, 



z ■ 



The reconstruction is carried out as described in section 2.3. This process requires solving 
the system ( [l6|) on each of the radial grid. If (3 is not too large, then (|l^) is an elliptic system 
of partial differential equations on the sphere 5^, mapping surjectively to the space of spin 2 
fields. We solve (|16|) by first substituting 

7 = S-'r, (74) 

where is defined spectrally by 

^''Yl = - + 2)(Z - l)r'^'Yil, I > 2, 

so (p^) becomes 
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run_160: 




Figure 8: Evolution of r/3 for < z < 55, in radial in-falling coordinates. Observe that the 
in-falling grid tracks the dynamical evolution. This simulation has timestep Az = 0.05, n = 
is the past horizon r = 2M and n = 256 represents future null infinity T"*". 



Note that the choice (^4^ gauges the I = 1 spherical harmonic components of 7 to zero — a 
similar but more expensive construction may be used if nonzero 7^=1 components are desired. 
The advantage of ( [75|) over (|l6|) is that the operator ICp in (|75| ) is close to the identity for small 
B := (5(3/(2 — div/3). The corresponding discretized problem is therefore well suited to iterative 
matrix methods. 



We use the conjugate gradient (CG) method |33|, which is an iterative method applicable to 
matrix problems of the form Ax = b with A symmetric positive definite. Accordingly we actually 
solve an associated self-adjoint equation obtained by applying to (^) the operator adjoint to 
that in (^) with respect to the norm on S"^. 

In order to compute the adjoint operator /Cj, notice first that £^ and /C^ are rea/-linear 
but not complex-linear, so the adjoint must be computed with respect to the real form of the 
inner product (|4^). Expanding T = X^i>2 T^'^Y';^, 4> = J2i>i4'''"^^im have the spectral 
representations 

g-^r = - -!)(; + 2)] -'/V'™y,L (76) 

l>2,m 

div</) = Yl + + (77) 

l>l,m 

and thus the real adjoints (9~^)^, div"^ are 

g-i^^ = - E m+m-i)r'^"<t>^"'Yi, (78) 

/>2,m 



div^(/ + i5) = 2Y,m + W f"'YL (79) 

l>l,m 

where /, g are real-valued functions. Consequently we may represent the adjoint /CT spectrally 
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by 



(/C^<^)^- = - J^Jl±^^{B<j> + B^)'-, (80) 



where B = 9/3/(2 - div/3) and / > 2. 

The equation ( p^ transformed into ( [75| ) gives an invertible equation 

ApF := /Cj/C/jF = /Cj(JS - K). (81) 

The right hand side of (|8l| ) may be computed exphcitly using the spectral representation and 
(|80|). The operator is symmetric and positive definite and close to the identity when con- 
sidered in the spectral representation. Thus the conjugate gradient algorithm may be applied 



to (81) and we see that the parameterisation of (|l6|) in terms of F (74) amounts to a precon- 
ditioner. Note CG requires not that the matrix of A/3 be given explicitly, but only that A/3^ 
can be evaluated for any spin 2 field We carry out this evaluation by a sequential process 



which computes the actions of ICp and /CT using the spectral representations of 9 and 



div, div^ (which are simple diagonal operators) combined with transformations to the point rep- 
resentation to evaluate multiplication terms like followed by projections back to the spectral 
representation. 

This scheme requires several transformations between representations of fields by their spin- 
weighted spherical harmonic coefficients and by their values on the S'^ grid. For example, the 
operator is a trivial multiplicative operator on spectral coefficients, whereas the products in 
the source terms are best calculated in the grid representation. Although evaluating the action 
of A/3 is thus numerically expensive, the expense is more than compensated for by the rapid 
convergence of the CG algorithm with this spectral preconditioning. 

The spectral representation has the further advantages that the solution is represented fewer 
unknowns (2(L + 1)^ — 8 compared to 4(L + 1)^ for the grid value representation), and gauge 
conditions which specify the / = 1 components of 7 (eg. ji=i = 0) can be directly implemented. 
It is possible to adapt the algorithm to allow for other NQS gauges (eg. Pi=i = 0), but this is 
numerically more expensive since (|^) must be solved 4 times at each sphere rather than once. 
Although such gauges have some geometric advantages their numerical implementation has 
not yet been considered. 

Using CG to solve for the spin 2 spherical harmonic coefficients of F turns out to be quite 
efficient, typically requiring fewer than 10 iterations for an S*^ grid of size N/2 x N = 16 x 32. 
On this size grid we resolve all components of F up to angular momentum L = N/2 — 1 = 15, 
so in this case we are solving for 2((L + 1)^ — 4) = 504 spectral coefficients. The scheme's 
effectiveness is due in part to having a good initial guess for F to use as the starting point of the 
CG iterations, namely the solution found for F on the 2-sphere at the previous radial position. 

The CG iterations finish when the error, measured by the sum of squares of spherical har- 
monic coefficients of the difference of the two sides of (|75|), is 10~^ times the size of the aliasing 
error in the source term. This aliasing error is the difference between the raw field values of the 
source term (which is necessarily calculated in the field value representation because it involves 
products and quotients) and its field values after projection into the subspace spanned by spin 2 
spherical harmonics. It provides an estimate of the error in the source term, and hence (because 
the operator IC/3 is close to the identity) it is reasonable to accept a solution of comparable 
accuracy. 

To ensure termination of the CG algorithm, other stopping criteria are also checked, but 
the relative error test is the usual termination cause and is found to work well in practice. For 
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example, it can result in a 10-fold improvement over letting the CG iterations run until the 
solution is determined to machine precision. 



4.5 Evolution 

Given /? on a null hypersurface Nz, we construct the time derivative dfj/dz by solving the 
hypersurface equations with seed /3, determining 7, v as outlined in the previous section, and 
then using formula (|l^) to evaluate d(3/dz. Let us write the result of this process as 

^ = B{(3,Uo) (82) 

where the operator B is determined by the value of /? on the hypersurface M and the initial 
conditions Uq = {Hq, Qq, Jq, Kq) at r = tq for the hypersurface equations. 

The evolution formula (|8^ ) provides the basis of the spacetime evolution algorithm, which 
simply incorporates ( p2[ ) into a standard 4th order Runge-Kutta algorithm. This approach is 
just the method of lines, treating the evolution equations as a very large system of ordinary 
differential equations for the (point) representation of the entire field P{z) = (3\j\f^- 

The method of lines, applied blindly in this manner, is generally prone to instabilities. Tests 
suggest the relative stability of the NQS code derives from the smoothing effects (a) of the 
convolution spline, and (b) of the spectral projection. The filtering implicit in the convolution 
spline is applied to P during the radial integration of the hypersurface equations, at each of the 
4 stages of the RK4 algorithm. It is not possible to turn off this radial filtering because the 
convolution splines for /3 are an essential part of the algorithm for evaluating the right hand side 

Smoothing of /5 in the angular directions is done explicitly, by projecting /3 onto the spin 1 
subspace with maximum angular momentum L or 2L/3 (the Orszag 2/3 rule, to eliminate 
quadratic aliasing) . This angular filtering is done after each of the 4 stages of the RK4 algorithm. 
Removal of the angular filtering results in very rapid disintegration of the evolution, which then 
typically lasts only a few RK4 steps. 

For simplicity, the 4 RK4 stages evolve /3 in the z direction in the {z, r) coordinates, along 
r = const. At the end of each full RK4 time step the key field /3 is interpolated onto the new 
radial grid (^), using a convolution spline based on values of P on the old grid. 

The RK4 time integration of P evolves field values on the (??, ip)-gnd. Equivalently, we could 
have evolved its spherical harmonic coefficients, of which there are only half as many. However, 
the amount of computation saved by doing so is insignificant in comparison to that required to 
evaluate dp/dz, so this choice is made for convenience. 

Likewise, the RK8 radial integration of the system of hypersurface equations uses the field 
value representation. In this case, however, it is found that projecting the fields onto their 
appropriate spherical harmonic subspaces during the integration is not required for either sta- 
bility or accuracy. There is a definite computational advantage in staying within the field value 
representation, since several relatively expensive 0{L^) projections are avoided. 

The first radial derivative of 7 is needed to evaluate dp/dz (p^. Grid values of dj/dr 
are calculated numerically as derivatives of convolution splines (in the radial direction) for the 



(f?, ip)-gnd values of 7, making use of formula (73) and the chain rule for derivatives. The radial 



derivative term P^logit appears in the hypersurface equations (^) and (|3^). Using equation 
( pOD and definition (|l5|), this term can be written as an expression involving only the 1st radial 
and angular derivatives of p. 

The program is normally run until the solution ceases to be well behaved. Blowup is detected 
by monitoring 2— div/3, which must remain everywhere positive. For the initial data that we have 
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used, the blowup has always occurred in the I = 2 modes of /?, at low values of n corresponding 
to r — 2M (see Figures Although the precise cause of the blowup is not yet understood, 
it is not a numerical instability, since it is unaffected by changes in the radial or timestep 
resolutions, nor does it appear to be primarily geometric, since most curvature scalars remain 
bounded. This suggests the blowup is a coordinate effect, probably arising from proximity to 
the future event horizon. 

For smooth initial data of intermediate strength (run_160), the evolution extends to z ~ 55. 
The final time varies with the strength of the initial data — see Table |3[ The evolution of an 
intermediate strength solution is shown in Figure ^, which plots the mean square or L?'{S'^) size 
of (3 at each radial sphere, for time < z < 55. Termination is caused by the blowup feature at 
low radius, which grows steadily from time z = 40 onwards. 



5 Accuracy tests 

The complexity of the NQS Einstein equations and the variety of algorithms employed in the 
code, make it problematic to prove rigorously that the numerical simulation accurately models 
the physics and geometry of the spacetime. Instead we rely on a range of tests to justify the 
reliability of the code, probing the numerical accuracy of the solutions through their convergence 
and geometric consistency. 

We consider here tests based on the numerical convergence of the solutions as algorithmic 
parameters are varied; and on the algebraic consistency of the numerical solutions. The con- 
sistency tests measure the constraint identities and the Trautman-Bondi mass decay formula 



p3| , 28]. Work in progress considers other tests, including comparisons with linearised theory, 
and with better known solutions such as Robinson- Trautman, Schwarzschild and Minkowski 
spacetimes in twisted NQS coordinates 0]. 

The resolution of the simulations is determined by three parameters: the spherical harmonic 
spectral limit L (or effective limit /max); the number of radial zones noo; and the time step Az. 
We shall examine in turn how the accuracy of a solution depends on each of these parameters. 

It is clear that numerical convergence can be estimated from the convergence properties of 
the key field /?. However, convergence of (3 guarantees only that the (limit) solution satisfies 
some system of equations, which may not coincide with the desired vacuum Einstein equations. 
(For example, the Einstein equations may have been incorrectly implemented). Thus, to assert 
that the correct equations have been solved, it is essential to provide independent tests of the 
correctness of the code. 

The most natural independent test is to compare the numerical solution with an explicitly 
known solution. Unfortunately the Schwarzschild metric @j is trivial in the NQS gauge and does 
not provide a useful comparison test, whilst the twisted shear- free metrics Q require boundary 
conditions which are more general than those available in the present version of the code. 

Instead we consider here another class of independent tests based on constraint relations. 
Such relations are typical of geometric equations arising in geometry and physics, which admit 
gauge and coordinate freedoms. Thus, we check the geometric consistency of the solution by 



evaluating r^G„„ and r^Gnm, using the constraint relations (26) and (p7|). Neither of these 
relations is used in generating the numerical solutions, and in theory these components should 
evaluate to zero. In practice, since each is a sum of terms having magnitude approximately 
1/3 1 ~ 1, the extent to which r^G„„, r'^Gnm evaluate to zero serves both to confirm the consistency 
of the numerical solution with the vacuum Einstein equations, and also to assess the accuracy 
of the solution. 
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The Trautman-Bondi mass decay formula provides another such test of geometric consistency, 
and of the accuracy of the solution near r = oo. This theoretical result leads to a relation between 
the asymptotic (r = oo) values and z-derivatives of the fields H, J, and K, and may be readily 
tested for our numerical solutions. 

In the following we discuss numerical solutions using three reference initial [5{z = 0) fields, 
which differ only by the scale factors given in Table ^. In each case the initial /3 consists of pure 
/ = 2, m = 2 spherical harmonics with equal strength odd and even parts, and radial profile 
given by a bump supported on 5 < r < 40. We use the terms weak, intermediate and strong to 
describe solutions generated using the three sizes of initial data. 

Another convenient measure of the strength of the gravitational field is the initial relative 
mass difference mB(0)/M— 1, between the initial Bondi mass of the numerical spacetime (cf. (|85|)) 
and the background Schwarzschild mass (M). Table I gives the initial relative mass differences 
for the three reference initial (5 fields. 



Table 3: Size of initial data sets 



Field strength: 


weak 


intermediate 


strong 


/?(0) scale factor: 


1 


4.48 


10 


mij(0)/M - 1: 


0.9472 X 10"^ 


0.1915 


0.9940 


Last z: 


61 


55 


51 



The qualitative conclusions of the error analysis of this section may be summarised as follows: 
The major determining factor in the overall accuracy is the spectral limit L. Truncating 
spherical harmonic coefficients beyond L has the effect of modifying the Einstein equations to 
a system for which the constraint identities are no longer valid, and thereby places a lower 
bound on the numerical accuracy. For L = 15 the weak and intermediate field solutions can be 
adequately resolved, but this is not sufficient to obtain adequate (beyond 10"'^) accuracy for the 
strong field simulation run_170. 

To suppress unstable quadratic aliasing effects, it is essential to use an Orszag 2/3 rule 
truncation. 

Within the bounds governed by the spectral limit L, accuracy can be improved by increasing 
the radial resolution Uoo- For the weak field solution, n^o = 1024 reduces the radial error 
contribution to the level of the spectral truncation error (see Figure ^4|(b)). 

For given resolutions L and Uoo, there is a range of values Az for which the simulation 
remains stable. Outside this range, the simulation follows the standard solution for some time, 
then rapidly blows up. The simulation is largely insensitive to the value of Az within the stable 
range, so Az may be chosen as large as possible, consistent with stable evolution. 



5.1 Dependence on spectral limit L 

Using our current hardware it is not generally feasible to run the code at L = 31, and L = 7 is 
too low to be of interest. The code is normally run at L = 15 resolution (giving a 16 x 32 (i9, ip) 
grid) with an anti-aliasing cutoff at Zmax = 10. 



Orszag ||l^, 41 1 observed that quadratic aliasing can be eliminated by periodically removing 



the upper 1/3 of the spectral bandwidth of a numerical solution. If fields contain only modes for 
which I < |L, then a quadratic product is band limited to / < |L. With a working bandwidth 
L, the modes for which L < I < |L become aliased onto the modes |L < I < L. Therefore, 
truncation at /max = 1-^ will remove quadratic aliasing contamination. 
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If no cutoff is used (ie. the full L = 15 resolution is retained) then the high /-modes of the 
intermediate strength simulations blow up at z ~ 8. The onset of instability (time until blow 
up) of the L = 15 simulations with no /max cutoff is largely independent of the time step and the 
radial resolution. This suggests that the effect of the nonlinear aliasing contamination is best 
regarded as changing the system of equations into a system which has unstable solutions. 

Because the nonlinear interactions in the NQS equations are predominantly quadratic, it is 
not surprising that the /max = 10 cutoff is sufficient for long term stability. The intermediate 
strength solution lasts until z = 55, when the code terminates for other reasons. 

Figure |9|(a) shows blow up of run_453, an L = 15 simulation of the intermediate field strength 
solution. The / = 15 modes show rapid growth beyond z = 6, indicating the instability of the 
aliasing feedback. Figure ^(b) shows the difference between run_453 and the stable simulation 
run_456, which has an /max = 10 cutoff. Until the onset of the high /-mode instability (ie. for 
z < 6) there is good agreement between the two simulations, with approximately lO"^'^ relative 
difference for / = 2 modes and 10~^ relative difference for / = 10 modes. 



run_453: |r|S|, z=0:S:l ruii_453 - run_456: |a(r|8)|j z=0:8:l 




Figure 9: Orszag's 2/3 rule is used to remove aliasing instability: (a) unstable evolution of 
the high /-modes of an L = 15 simulation (no anti-aliasing cutoff); (b) difference between an 
unstable L = 15 simulation (no cutoff) and a stable simulation with an /max = 10 cutoff. Each 
/-bin contains a radial plot (linear in n, with n = 0, . . . ,noo) of the square root of the sum of 
the squares of the (/, m)-components for fixed / with m = —I, ...,/. These simulations have 
noo = 512 and Az = 0.05. 



The spectral limit L is critical in determining the relation between gravitational field strength 
and simulation accuracy. This can be appreciated by observing the decay rate of the /-spectrum 



of P, as in Figure |10[ By extrapolation, the error introduced by the anti-aliasing cutoff at 
^max = 10 should be no more than the / = 10 coefficient, and a relative error estimate follows by 
comparing the / = 10 and the / = 2 coefficients. 

Figures |l^(a) and |l^(b) show the dramatic difference in decay rates of the /-modes of (3 for 
weak and strong fields. Assuming an /max = 10 cutoff, it is evident that the relative error is at 
most 10~^ for weak field simulations, and about 10~^ for strong field simulations. 

From the observed decay rate of the /-modes of /? for a given field strength, it is possible to 
estimate the resolution L required to achieve a prescribed accuracy. Thus although we cannot 
directly investigate the behaviour of errors with varying spectral limit L (due to hardware 
constraints), we can still investigate spectral resolution effects by altering the P field strength. 

Figures 11(a) and 0(b) show the effect of field strength (weak, intermediate, strong) on the 
constraint quantities r'^Gnn and r'^Gnm- The parameters for these simulations are /max = 10, 
rioo = 256 and Az = 0.05. The four curves in each band are times z = 10, 20, 30, 40. There is no 
significant z dependence of either Gnn or Gnm until within about 5M of the final blow up time. 
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ruii_150: | | z=10 



run_170: | | z=10 




123456789 10 123456789 10 

(a) I (b) ! 

Figure 10: Spectral resolution and field strength: (a) well resolved weak field with fast ^-mode 
decay; (b) poorly resolved strong field with slow decay of /-modes (see Table ^ for field strength 
details) . 



run_150, 160 & 170: I r^Gj^^ l,sio z=10:40:10 run_150, 160 & 170: I r^G^^^ | , z=10:40:10 




64 128 192 256 64 128 192 256 



(a) n (b) n 

Fi gure 11: Effect of spectral resolution on constraint quantities (a) Ir^Gnnl^a; (b) \r'^Gnm\s'^i 
at times z = 10,20,30,40, for strong (top 4 curves), intermediate (middle 4 curves) and weak 
(bottom 4 curves) fields. 



The second Bianchi identity implies the conservation law G^^^ = 0, which leads to a radial 
system of equations for Gnn,Gnm with sources linear in the hypersurface Einstein tensor com- 
ponents Gil, Gim, Gin, Gmm- Thus Gnn,Gnm give a measure of the accumulated error in the 
hypersurface equations in the radial direction. This provides some explanation of the structure 
of the Gnn, Gnm graphs, particularly for the strong field solution: the numerical solution of the 
hypersurface equations will have greatest error in the region where the fields are strongest, in this 
case the range 64 < n < 128, and this is precisely the region of greatest increase in Gnm- 



5.2 Dependence on radial grid resolution rioo 

The radial regridding and interpolation of /3, the radial differentiation of 7, and the RK inte- 
gration of the hypersurface equations are all formally 8th order accurate. Figure 12 shows that 
this is consistent with the observed convergence of /? on increasing the radial resolution. 

The constraint quantities Gnn and Gnm also exhibit some convergence effects. Figures p^(a) 
and 13(b) show show significant improvement between rioo = 256 and 512, but little between 
512 and 1024. The form of the 77-00 — 1024 curve indicates that errors at the highest radial 
resolution are dominated by errors associated with the spectral truncation. 
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(run_352 & 452) - run_552: | i5(r(3) | , z=l 



in 




64 12B 192 256 



Figure 12: Convergence of (3 with increasing radial resolution: weak field solutions with = 
256, 512 compared to noo = 1024. The error decreases by approximately a factor of 2^ on 
doubling the radial resolution. 



run_352, 452 & 552: I r^G J z=l ruii_352, 452 & 552: I r^G | z=l 




64 128 192 256 64 128 192 256 



(a) n (b) n 

Figure 13: Effect of radial resolution on weak field constraint quantities (a) \r'^Gnn\s^^ i^) 
\r'^Gnm\s^- III each case the three curves are noo = 256,512,1024 (top, middle and bottom 
curves respectively). 



5.3 Dependence on time step Az 

Although the solution algorithm is formally 4th order accurate in the time direction, at the 
typical resolutions at which the code is run, the RK4 errors are completely dominated by errors 
arising from the spectral truncation L and/or the radial discretisation noo- This is illustrated 
by Figure 0(a), which shows no significant difference in the constraint quantity G„„ between 
Az = 0.1 and Az = 0.05, with noo = 256. However, when the solution is better resolved in 



the radial direction, a small effect can be observed, cf. Figure 14(b), where noo = 1024. Figure 



15 compares r/3 for runs with Az = 0.1,0.05,0.025 and noo = 512 and again shows only minor 
improvements from decreasing Az. 

Consequently, Az is optimally chosen as large as possible, subject to resulting in stable 
evolution. For noo = 256 and L = 15 with an anti-aliasing luia.x = 10 cutoff, the evolution is 
stable for Az = 0.1 and unstable for Az = 0.2, which blows up at time z = 25, after 125 RK4 
steps. 



31 




64 128 192 256 128 256 384 512 640 768 896 1024 



(a) n (b) n 

Fi gure 14: Effect of time step resolution on the constraint quantity \r'^Gnn\s'^ ^'^^ the weak field 
solution: (a) Hoo = 256, Az = 0.1, 0.05: the error is dominated by the radial discretisation error 
for n < 192 and by the spectral truncation error for n > 192. Refining Az produces no ap- 
preciable improvement in the solution, (b) no© = 1024, Az = 0.1,0.05: the radial discretisation 
error is small enough that the RK4 integration error can be observed. For n > 700 the error is 
dominated by the spectral truncation, resulting in the same tail as in (a), while for n < 700 the 
constraint improves in places, consistent with a factor of 16 decrease in the error. 



{run_442 & 452) - run_462: \6(T|S)\,^^^ z=l 




128 256 384 512 



Figure 15: Convergence of f3 with decreasing time step: weak field solutions for Az = 0.1,0.05, 
compared against Az = 0.025. Where the error is not dominated by the radial discretisation 
error, the curves show a decrease in error which is consistent with 4th order convergence. 



5.4 Energy and aisymptotic decay tests 

The Hawking mass 



of a 2-surface S reduces in the NQS gauge to 

mniz, r) = \r{l-^ jf ^ HJ^ . (84) 

mH{z,r) provides an easily computed quantity representing the "quasi-local" mass contained 
within the sphere {z,r), and has asymptotic limit equal to the Bondi mass 

'niB{z) = lim mHir,z). (85) 
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The Bondi mass is easily computed numerically, by mB{z) = mnin = n^o^z). Figure p^(a) 
shows the Hawking mass plotted against the radial coordinate, for times z = 0, 1, . . . , 60. There 
are several features of interest in this plot: the limit Bondi mass (Figure |T6|(b)) decays in time, 
reflecting the Trautman-Bondi mass loss formula (^); the energy is radiated in bursts, reflecting 
near-linear behaviour dominated by pure I = 2 modes; the Hawking and Bondi masses decay to 
background black hole mass M = 1 at late times, suggesting that in this example, almost all the 
gravitational radiation has been scattered to and essentially none will be absorbed by the 
black hole; and finally, the rapidly growing feature about n = 20 at late times in Figure does 
not affect the Hawking mass. 

The Trautman-Bondi mass loss formula |^, 14, 28| 



az 



lim 

IGvr r— »oo 



H\K\ 



(86) 



provides another test of the geometric consistency of the solution, particularly near null infinity. 
By comparing the numerical derivative dm s/dz with the computed value of the right hand side 
(evaluated at n = rioo), we may construct the error ^m-s — RHS (p6|) . Figure |l^(b) plots this 



error against time z, suggesting that the asymptotic (r 
about 0.00001%. 



oo) fields of run_160 are accurate to 



run_160: Bondi mass 



run_160: Hawking mass z — 0:55:1 




(a) 

Figure 16: Mass functions: (a) Hawking mass for times 2: = 0, 1, ... , 55; (b) Bondi mass. 



run_160: Bondi mass decay rate 



run_160: Error in mass decay rate 
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Figure 17: The Trautman-Bondi mass loss formula as test of numerical accuracy at r = cxd: 
(a) Bondi mass decay rate; (b) error in the mass decay formula, given by LHS(^^ — RHS (^6|) . 



References 



33 



[1] S. Abarbanel, D. Gottlieb, and M. H. Carpenter. On the removal of boundary errors caused 
by Runge-Kutta integration of nonlinear partial differential equations. SI AM J. Sci. Comp., 
17:777-782, 1996. 

[2] R. Balean. The null-timelike boundary problem for the linear wave equation. PhD thesis, 
University of New England, 1996. 

[3] R. Balean. The null-timelike boundary problem for the linear wave equation. Commun. 
P.D.E., 22:1325-1360, 1997. 

[4] R. Bartnik. Quasi-spherical metrics and prescribed scalar curvature. J. Diff. Geom., 37:31- 
71, 1993. 

[5] R. Bartnik. Einstein equations in the null quasi-spherical gauge. Class. Quant. Gravity, 
14:2185-2194, 1997. 

[6] R. Bartnik. Einstein equations in the null quasi-spherical gauge: progress report. In 
D. Wiltshire, editor, ACGRGl Proceedings, Adelaide February 1996, pages 25-38, 1997. 
http : //www . physics . adelaide . edu . au/ASGRG/ ACGRGl/. 

[7] R. Bartnik. Shear-free null quasi-spherical spacetimes. J. Math. Phys., 38:5774-5791, 1997. 

[8] R. Bartnik. Interaction of gravitational waves with a black hole. In T. Bracken and D. De 
Wit, editors. Proceedings, Xllth Int'l. Congress Math. Phys., Brisbane, July 1997. Interna- 
tional Press, to appear 1999. 

[9] R. Bartnik and A. H. Norton. NQS data explorer. Technical report. University of Canberra, 
1997. 

http : / / gular . Canberra . edu . au/relat ivity . html. 

[10] R. Bartnik and A. H. Norton. Numerical solution of the Einstein equations. In A. Gill 
J. Noye, M. Teubner, editor. Computational Techniques and Applications: CTAC97, pages 
91-98, 1998. 

[11] N. Bishop, R. Gomez, L. Lehncr, M. Maharaj, and J. Winicour. High-powered gravitational 
news. Phys. Rev. D, 56:6298, 1997. 

[12] N. T. Bishop. Some aspects of the characcristic initial value problem in numerical relativity. 
In C. J. Clarke and R. d'Inverno, editors. Prospects in Numerical Relativity II, pages 20-33. 
Cambridge, 1992. 

[13] S. Bonazzola, E. Gourgoulhon, and J. -A. Marck. Spectral methods in general relativistic 
astrophysics, gr-qc/9811089, 1998. 

[14] H. Bondi. Gravitational waves in general relativity. Nature, 186:535, 1960. 

[15] J. P. Boyd. The choice of spectral functions on a sphere for boundary and eigenvalue 
problems: A comparison of Chebyshev, Fourier and asssociated Legendre expansions. Mon. 
Weather Rev., 106:1184-1191, 1978. 

[16] J. P. Boyd. Chebyshev and Fourier Spectral Methods. Lecture Notes in Engineering, vol 49. 
Springer, 1989. 



34 



[17] G. L. Browning, J. J. Hack, and P. N. Swarztrauber. A comparison of three numerical 
methods for solving differential equations on the sphere. Mon. Weather Rev., 117:1058- 
1075, 1989. 

[18] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid 
Mechanics. Series in Computational Physics. Springer, 1988. 

[19] S. Chandrasekhar. The Mathematical Theory of Black Holes. Oxford, 1984. 

[20] R. D'Inverno. Introducing Einstein's Relativity. Oxford, 1992. 

[21] M. Eastwood and P. Tod. Edth - a differential operator on the sphere. Math. Proc. Camb. 
Phil. Soc, 92:317-330, 1982. 

[22] J.A.H. Futterman, F.A. Handler, and R.A. Matzner. Scattering from black holes. Cam- 
bridge, 1988. 

[23] J. N. Goldberg, A. J. Macfarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan. 
Spin-s spherical harmonics and S. J. Math. Phys., 8:2155-2161, 1967. 

[24] R. Gomez, P. Papadopoulos, and J. Winicour. Null cone evolution of axi-symmetric vacuum 
space-times. J. Math. Phys., 35:4184-4203, 1994. 

[25] R. Gomez and J. Winicour. Asymptotics of gravitational collapse of scalar waves. J. Math. 
Phys., 33:1445, 1992. 

[26] R. Gomez, J. Winicour, and R. Isaacson. Evolution of scalar fields from characteristic data. 
J. Comp. Phys., 98:11, 1992. 

[27] R. Gomez, J. Winicour, and P. Papadopoulos. The eth formalism in numerical relativity. 
Class. Quantum Grav., 14:977-990, 1997. 

[28] M. G. J. van der Berg H. Bondi and A. W. K. Mctzner. Gravitational waves in general 
relativity VH. Proc. Roy. Soc. Lond., A269:21-51, 1962. 

[29] R. W. Hamming. Digital Filters. Prentice-Hall, 1977. 

[30] S. W. Hawking and G. R. Ellis. The large-scale structure of spacetime. Cambridge UP, 
Cambridge, 1973. 

[31] A. C. Hearn. REDUCE Users Manual, Version 3.6. RAND Publication CP78 (rev. 9/96), 
1996. 

[32] R. Jakob-Chien and B. K. Alpert. A fast spherical filter with uniform resolution. J. Comp. 
Phys., 136:580-584, 1997. 

[33] C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element 
Method. Cambridge, 1987. 

[34] P. E. Merilees. The pseudo-spectral approximation applied to the shallow water equations 
on a sphere. Atmosphere, 11:13-20, 1973. 

[35] V. Moncrief. Gravitational perturbations of spherically symmetric systems. Ann. Phys., 
88:323-342, 1974. 



35 



E. Newman and R. Penrose. An approach to gravitational radiation by a method of spin 
coefficients. J. Math. Phys., 3:566-578, 1962. 

E. T. Newman and T. Unti. A class of null flat-spce coordinate systems. J. Math. Phys., 
4:1467-1469, 1963. 

A. H. Norton. Finite difference operators for pdcs, based on sampling kernels for spline 
quasi-intcrpolation. Technical report. School of Mathematics, UNSW, 1992. 

A. H. Norton. Spectral collocation methods for solution of Einstein's equations in null 
quasi-spherical coordinates. In D. L. Wiltshire, editor. Proceedings of the First Aus- 
tralasian Conference on General Relativity, Adelaide February 1996, pages 39-47, 1997. 
http : //www . physics . adelaide . edu . au/ASGRG/ACGRGl/. 

A. H. Norton. The use of convolution splines in the NQS Einstein evolution code. In A. Gill 
J. Noye, M. Teubner, editor. Computational Techniques and Applications: CTAC97, pages 
473-480, 1998. 

S. A. Orszag. On the elimination of aliasing in finite difference schemes by filtering high 
wave number components. J. Atmosph. Sci., 28:1074, 1971. 

S. A. Orszag. Fourier series on spheres. Monthly Weather Rev, 102:56-75, 1974. 

R. Penrose. The geometry of impulsive gravitational waves. In L. O'Raifaertaigh, editor. 
General Relativity: in honour of J.L.Synge. Oxford, 1972. 

R. Penrose and W. Rindler. Spinors and Space-time, vol. I, II. Cambridge UP, 1984,1986. 

D. A. Prager and A. W. C. Lun. Computational methods in the physical interpretation of 
Robinson-Trautman spacetimes. In D. L. Wiltshire, editor. First Australasian Conference 
on General Relativity and Gravitation, 1996. 

P. J. Prince and J. R. Dormand. High order embedded Runge-Kutta formulae. J. Comp. 
Appl. Math., 7:67-75, 1981. 

T. Regge and J. A. Wheeler. Stability of a Schwarzschild singularity. Phys. Rev., 108:1063- 
1069, 1957. 

A. D. Rendall. Reduction of the characteristic initial value problem to the Cauchy problem 
and its application to the Einstein equations. Proc. Roy. Soc. Lond., A427:221-239, 1990 
MPI preprint 438, 1989. 

I. Robinson and A. Trautman. Spherical gravitational waves in general relativity. 
Proc. Roy. Soc. Lond. A, 270:103-126, 1962. 

R. K. Sachs. Gravitational waves in general relativity, VIII. waves in asymptotically fiat 
space-time. Proc. Roy. Soc. Lond. A, A270:103-126, 1962. 

R. K. Sachs. Gravitational waves in general relativity VIII. waves in asymptotically fiat 
space-time. Proc. Roy. Soc. Lond., A270: 103-126, 1962. 

L. J. Schocnberg. Cardinal Spline Interpolation, volume 12 of Regional Conference Series 
in Applied Mathematics. SIAM, 1973. 



36 



[53] D. Singleton. Robinson- Trautman solution of Einstein's equations. PhD thesis, Monash 
University, 1990. 

[54] L. Smarr and Jr. J. W. York. Kinematical conditions in the construction of spacetime. 
Phys. Rev. D, 17:2529-2551, 1978. 

[55] G. Starius. Composite mesh difference methods for elhptic and boundary value problems. 
Numer. Math., 28:243-248, 1977. 

[56] K. Stellmacher. Ausbretungsgesetze fiir charakteristische Singularitaten der Gravitations- 
gleichungen. Math. Annalen, 145:740-783, 1938. 

[57] F. Stenger. Numerical Methods Based on Sine and Analytic Functions. Springer, 1993. 

[58] P. N. SwarztraTibcr. On the spectral approximation of discrete scalar and vector functions 
on the sphere. SIAM J. Numer. Anal, 16:934-949, 1979. 

[59] P. N. Swarztraubcr. The approximation of vector functions and their derivatives on the 
sphere. SIAM J. Numer. Anal., 18:191-210, 1981. 

[60] P. N. Swarztrauber. Software for the spectral analysis of scalar and vector functions on the 
sphere. In Large Scale Scientific Computation. Academic Press, 1984. 

[61] P. N. Swarztraubcr. Spectral transform methods for solving the shallow-water equations 
on the sphere. Mon. Weather Rev,, 124:730, 1996. 

[62] A. Trautman. King's College lecture notes on general relativity. May-June 1958. 

[63] A. Trautman. Radiation and boundary conditions in the theory of gravitation. Bull. Acad. 
Pol. Sci., Serie sci. math., VI:407-412, 1958. 

[64] A. Trautman. Conservation laws in general relativity. In L. Witten, editor. Gravitation: 
an introduction to current research. Wiley, 1962. 

[65] F. J. Zerilli. Effective potential for even parity regge-wheeler gravitational perturbation 
equations. Phys. Rev. Lett., 24:737-738, 1970. 

[66] H. Miiller zum Hagen. Characteristic initial value problem for hyperbolic systems of second 
order differential equations. Ann. Inst. Henri Poincare, 53:159-216, 1990. 



37 



